## BELL (2024) ESTIMATOR FOR LISS

# Create Dummies
clean_LISS_data_Bell <- fastDummies::dummy_cols(clean_LISS_data, select_columns = c("wage_q"))

# Group Controls
wage_controls <- colnames(clean_LISS_data_Bell)[grep("^wage_q_", colnames(clean_LISS_data_Bell))]
wage_controls <- wage_controls[!wage_controls %in% "wage_q_1"] # omit lowest wage category as base

controls <- paste(c(wage_controls), collapse = " + ")
regression_equation_LISS <- paste("education ~ 0 + high_meaning + high_flex + high_telecommute +", controls)

# First Stage
first_stage <- lm(formula = regression_equation_LISS, data = clean_LISS_data_Bell)
rss_full <- sum(residuals(first_stage)^2)
n = nobs(first_stage)
k = qr(model.matrix(first_stage))$rank

# Save first-stage for Table
first_stage_liss <- first_stage

# Compute Coefficients + Anderson-Rubin bounds
pi <- coef(first_stage)["high_meaning"]
delta <- coef(first_stage)
cov_matrix <- vcov(first_stage)

c_pipi <- cov_matrix["high_meaning", "high_meaning"]
c_deltadelta <- diag(cov_matrix)
c_delpi <-  cov_matrix["high_meaning", ]
crit <- 1.96

a <- pi^2 - (crit^2 * c_pipi)
b <- 2 * crit^2 * c_delpi - 2*delta*pi
c <- ( delta^2 - (crit^2)*c_deltadelta)

lb <- (- (-b + sqrt(b^2 - 4 * a * c)) / (2 * a))
ub <- (- (-b - sqrt(b^2 - 4 * a * c)) / (2 * a))
beta <- -delta/pi

# Nested First Stage
regression_equation_nest <- paste("education ~ 0 + high_flex + high_telecommute +", controls)
first_stage_nest <- lm(formula = regression_equation_nest, data = clean_LISS_data_Bell)
rss_nest <- sum(residuals(first_stage_nest)^2)

# Compute partial F 
partial_F = (rss_nest - rss_full) / (rss_full / (n - k))

# Linear Regressions (base and with productivity controls) + Bell Estimates
base_LISS <- lm(high_meaning ~ high_flex + high_telecommute, data = clean_LISS_data)
prod_LISS <- lm(high_meaning ~ high_flex + high_telecommute + education, data = clean_LISS_data)

# Confidence intervals
conf_base_LISS <- confint(base_LISS, level = 0.95)
conf_prod_LISS <- confint(prod_LISS, level = 0.95)

# Create Table
base_flex_LISS <- data.frame(conf_base_LISS["high_flex",1], base_LISS$coefficients["high_flex"], conf_base_LISS["high_flex",2]) 
base_flex_LISS[] <- lapply(base_flex_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})

base_tel_LISS <- data.frame(conf_base_LISS["high_telecommute",1], base_LISS$coefficients["high_telecommute"], conf_base_LISS["high_telecommute",2])
base_tel_LISS[] <- lapply(base_tel_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})

prod_flex_LISS <- data.frame(conf_prod_LISS["high_flex",1], prod_LISS$coefficients["high_flex"], conf_prod_LISS["high_flex",2]) 
prod_flex_LISS[] <- lapply(prod_flex_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})

prod_tel_LISS <- data.frame(conf_prod_LISS["high_telecommute",1], prod_LISS$coefficients["high_telecommute"], conf_prod_LISS["high_telecommute",2])
prod_tel_LISS[] <- lapply(prod_tel_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})

Bell_flex_LISS <- data.frame(lb["high_flex"], beta["high_flex"], ub["high_flex"])
Bell_flex_LISS[] <- lapply(Bell_flex_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})

Bell_tel_LISS <- data.frame(lb["high_telecommute"], beta["high_telecommute"], ub["high_telecommute"])
Bell_tel_LISS[] <- lapply(Bell_tel_LISS, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})


## Open a connection to a file
file_conn <- file(file.path(graph_dir, "Proxy_LISS.txt"), "w")

# Write the LaTeX code to the file
writeLines("\\begin{tabular}{llll}", file_conn)
writeLines("\\hline \\hline", file_conn)
writeLines("              & Base        & Productivity Controls & Bell Proxy  \\\\", file_conn)
writeLines("\\hline", file_conn)

# Write point estimates for TELECOMMUTE
writeLines(paste("Telecommuting", " & ",
                 base_tel_LISS[1,2], " & ", 
                 prod_tel_LISS[1,2], " & ", 
                 Bell_tel_LISS[1,2], " \\\\"), file_conn)

# Write confidence intervals for TELECOMMUTE
writeLines(paste("   \\textit{Conf. Int.}           & ", 
                 "(", base_tel_LISS[1,1], ",", base_tel_LISS[1,3], ")", " & ",
                 "(", prod_tel_LISS[1,1], ",", prod_tel_LISS[1,3], ")", " & ",
                 "(", Bell_tel_LISS[1,1], ",", Bell_tel_LISS[1,3], ")", " \\\\[0.5em]"), file_conn)

# Write point estimates for FLEXIBILITY
writeLines(paste("High Flexibility", " & ",
                 base_flex_LISS[1,2], " & ", 
                 prod_flex_LISS[1,2], " & ", 
                 Bell_flex_LISS[1,2], " \\\\"), file_conn)

# Write confidence intervals for FLEXIBILITY
writeLines(paste("   \\textit{Conf. Int.}           & ", 
                 "(", base_flex_LISS[1,1], ",", base_flex_LISS[1,3], ")", " & ",
                 "(", prod_flex_LISS[1,1], "," , prod_flex_LISS[1,3], ")", " & ",
                 "(", Bell_flex_LISS[1,1], ",", Bell_flex_LISS[1,3], ")", " \\\\[0.25em]"), file_conn)

# Write partial F row
writeLines("\\hline", file_conn)
writeLines(paste("Partial F:    &             &                       &", round(partial_F, 3),          "\\\\"), file_conn)
writeLines("\\hline", file_conn)
writeLines("\\hline", file_conn)
writeLines("\\\\", file_conn)

# End table and close file connection
writeLines("\\end{tabular}", file_conn)
close(file_conn)

rm("a", "b", "c_pipi", "c_deltadelta", "c_delpi", "crit", "cov_matrix", "pi", "delta", "n",  "k", "rss_full", "controls")
